function [gbestx,gbestfitness,gbesthistory] = GEO (N,nvars,ub,lb,~,~,maxiter,fun,FuncId)
% Golden Eagle Optimizer (GEO)
%% initialization
options.PopulationSize = N;
options.MaxIterations  = maxiter;
options.AttackPropensity = [0.5 ,   2];
options.CruisePropensity = [1   , 0.5];
PopulationSize = options.PopulationSize;
MaxIterations = options.MaxIterations;
FEs=0;
MaxFEs=1e4*nvars;
gbestfitness=inf;
x = lb + rand (PopulationSize,nvars) .* (ub-lb);
FitnessScores = fun (x',FuncId);
[fmin,ind]=min(FitnessScores);
if fmin<=gbestfitness
    gbestfitness=fmin;
    gbestx=x(ind,:);
end
FEs=FEs+PopulationSize;
gbesthistory(1:FEs)=gbestfitness;
% solver-specific initialization
FlockMemoryF = FitnessScores;
FlockMemoryX = x;

AttackPropensity = linspace (options.AttackPropensity(1), options.AttackPropensity(2), MaxIterations);
CruisePropensity = linspace (options.CruisePropensity(1), options.CruisePropensity(2), MaxIterations);

%% main loop

for CurrentIteration = 1 : MaxIterations
    
    % prey selection (one-to-one mapping)
    DestinationEagle = randperm (PopulationSize)';
    
    % calculate AttackVectorInitial (Eq. 1 in paper)
    AttackVectorInitial = FlockMemoryX (DestinationEagle,:) - x;
    
    % calculate Radius
    Radius = VecNorm (AttackVectorInitial, 2, 2);
    
    % determine converged and unconverged eagles
    ConvergedEagles = sum (Radius,2) == 0;
    UnconvergedEagles = ~ ConvergedEagles;
    
    % initialize CruiseVectorInitial
    CruiseVectorInitial = 2 .* rand (PopulationSize, nvars) - 1; % [-1,1]
    
    % correct vectors for converged eagles
    AttackVectorInitial (ConvergedEagles, :) = 0;
    CruiseVectorInitial (ConvergedEagles, :) = 0;
    
    % determine constrained and free variables
    for i1 = 1 : PopulationSize
        if UnconvergedEagles (i1)
            vConstrained = false ([1, nvars]); % mask
            idx = datasample (find(AttackVectorInitial(i1,:)), 1, 2);
            vConstrained (idx) = 1;
            vFree = ~vConstrained;
            CruiseVectorInitial (i1,idx) = - sum(AttackVectorInitial(i1,vFree).*CruiseVectorInitial(i1,vFree),2) ./ (AttackVectorInitial(i1,vConstrained)); % (Eq. 4 in paper)
        end
    end
    
    % calculate unit vectors
    AttackVectorUnit = AttackVectorInitial ./ VecNorm (AttackVectorInitial, 2, 2);
    CruiseVectorUnit = CruiseVectorInitial ./ VecNorm (CruiseVectorInitial, 2, 2);
    
    % correct vectors for converged eagles
    AttackVectorUnit(ConvergedEagles,:) = 0;
    CruiseVectorUnit(ConvergedEagles,:) = 0;
    
    % calculate movement vectors
    AttackVector = rand (PopulationSize, 1) .* AttackPropensity(CurrentIteration) .* Radius .* AttackVectorUnit; % (first term of Eq. 6 in paper)
    CruiseVector = rand (PopulationSize, 1) .* CruisePropensity(CurrentIteration) .* Radius .* CruiseVectorUnit; % (second term of Eq. 6 in paper)
    StepVector = AttackVector + CruiseVector;
    
    % calculate new x
    x = x + StepVector;
    
x=min(x,ub);
x=max(x,lb);

    
    % calculate fitness
    FitnessScores = fun (x',FuncId);
    [fmin,ind]=min(FitnessScores);
    if fmin<=gbestfitness
        gbestfitness=fmin;
        gbestx=x(ind,:);
    end
    FEs=FEs+PopulationSize;
    gbesthistory(FEs-PopulationSize+1:FEs)=gbestfitness;
    
    % update memory
    UpdateMask = FitnessScores < FlockMemoryF;
    FlockMemoryF (UpdateMask) = FitnessScores (UpdateMask);
    FlockMemoryX (UpdateMask,:) = x (UpdateMask,:);
    if FEs>=MaxFEs
        break;
    end
    fprintf("GEO 第%d次评价，最佳适应度 = %e\n",FEs,gbestfitness);
end
if FEs<MaxFEs
    gbesthistory(FEs+1:MaxFEs)=gbestfitness;
else
    if FEs>MaxFEs
        gbesthistory(MaxFEs+1:end)=[];
    end
end
end
function N = VecNorm (A, p, dim)  
N = nthroot(sum(A.^p,dim),p); 
end

